#coding=utf-8

import cv2
import numpy as np
import math
import matplotlib.pyplot as plt
from numpy import unravel_index

MIN_VAL = 20.0                     # 该值用于拉伸；
MI = 2
limitVal = 50.0

def CalcCSF(v):
    """
    计算CSF;
    :param v:
    :return:
    """
    # a = 75
    # b = 0.2
    # c = 0.8
    # K = 34.05
    a = 1.0
    b = 0.2
    c = 0.8
    K = 1.0
    # K = 30
    return 1.0 * a * math.pow(v, c) * math.exp(-1.0 * b * v)/K

def CalcDetailReserve(mtf, f):
    """
    计算细节保留值；
    :param spatialImg:
    :return:
    """
    nLen = len(f)
    csf = np.array([CalcCSF(v) for v in f])
    dr = np.sum(mtf * csf)/np.sum(csf)
    # dr =(dr - MIN_VAL) * 1.0/(100.0 - MIN_VAL)                # 对比度拉伸；
    # if dr * 100.0 < limitVal:
    #     dr = math.pow(1.0 * dr * 100.0/limitVal, MI) * limitVal
    # else:
    #     dr = limitVal + math.pow(1.0 * (100.0 * dr - limitVal)/(100.0 - limitVal), 1.0/MI) * (100.0 - limitVal)
    return dr

def CalcLinearSpatialAndMtf(mtfImg):
    height, width = mtfImg.shape
    row = int(0.5 * (height + 1))
    col = int(0.5 * (width + 1))
    nbins = row
    amp = np.zeros((1, nbins), dtype="float64")
    stdamp = np.zeros((1, nbins), dtype="float64")
    freq = np.zeros((row, col), dtype = "int")
    maxVal = math.sqrt(row * row + col * col)
    for r in range(1, row + 1):
        for c in range(1, col + 1):
            freq[r-1][c-1] = int(1.0 * math.sqrt(r * r + c * c)/maxVal * nbins + 0.5) -1
    for r in range(1, row + 1):
        for c in range(1, col + 1):
            ind = freq[r - 1][c - 1]
            amp[0][ind] = amp[0][ind] + mtfImg[r - 1][c - 1]
            stdamp[0][ind] = stdamp[0][ind] + 1
    # for r in range(nbins):
    #     amp[0][r] = 1.0 * amp[0][r]/stdamp[0][r]
    f = 1.0 * np.array(range(nbins))/nbins * math.sqrt(0.5)
    return amp, stdamp, f

def CalcSpatialImg(img):                    # 空间频率;
    """
    计算频谱图；
    :param path:
    :return:
    """
    # img = cv2.imread("test.jpg", 0)
    dftImg = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dftImg_shift = np.fft.fftshift(dftImg)
    # dftImg_shift = np.array(dftImg_shift, dtype="uint8")
    result = 20 * np.log(cv2.magnitude(dftImg_shift[:, :, 0], dftImg_shift[:, :, 1]))
    result = result / np.max(result)
    return result

if __name__ == "__main__":
    img = cv2.imread("test.jpg", 0)
    img = cv2.medianBlur(img, 5)
    dftImg = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    dftImg_shift = np.fft.fftshift(dftImg)
    # dftImg_shift = np.array(dftImg_shift, dtype="uint8")
    result = 20 * np.log(cv2.magnitude(dftImg_shift[:, :, 0], dftImg_shift[:, :, 1]))
    result = result/np.max(result)
    mtf,stdMtf, f = CalcLinearSpatialAndMtf(result)
    dr_1 = CalcDetailReserve(mtf, f)
    dr_2 = CalcDetailReserve(stdMtf, f)
    dr = dr_1/dr_2 * 100.0
    startDr = 30.0
    endDr = 58.0
    dr = startDr + 100.0/(endDr - startDr) * (dr - startDr)
    print(dr)

    plt.subplot(122)
    plt.imshow(result, cmap ='gray')
    plt.title('result')
    plt.axis('off')
    plt.show()